Massively Parallel Simulations of Chain Formation and Restructuring Dynamics in a Magnetorheological Fluid

Magnetorheological fluids consist of micron sized iron particles mixed in a carrier fluid, and are commonly used in adaptive dampers. Current simulations of MR fluids have been limited to thousands of particle and have been unable to simulate a a practical fluid volume (∼mm3) with a high solids loading (∼ 25 vol %). In this paper, we use NVIDIA’s CUDA programming environment to simulate over one million particles. Using these simulations, we can dynamically simulate chain formation and restructuring in a practical millimeter scale fluid volume with realistic solids loading. The chain structures can be characterized in terms of extent and number, as well as obtaining physical metrics such as yield stress. Nomenclature χ Magnetic susceptibility ∆t Integration time step γ̇ Shear rate η Fluid viscosity r̂ Position unit vector ai Acceleration of particle i ex Unit vector in x direction Fi j Contact forces between particles i and j Fi Hydrodynamic force on particle i Fi j Magnetic forces between particles i and j H0 External applied field M Magnetization m Magnetic moment r Position vector v∞ Fluid velocity vi Velocity of particle i μ Magnetic permeability φ Volume fraction τ Shear stress τy Yield stress a Particle radius c Connectivity c/c0 Non-dimensional connectivity c0 Characteristic connectivity value Cd Coefficient of drag Cp Coefficient of polarization h Height of volume element k Repulsion constant L Mean chain length L/L0 Non-dimensional chain length L0 Characteristic chain length mi Mass of particle i N Number of particles n Spring power Qi j Repulsion strength r Magnitude of position vector Rc Cutoff radius Re Reynolds Number 1 Copyright c © 2011 by ASME INTRODUCTION Magnetorheological (MR) fluids and electrorheological (ER) fluids are non-Newtonian fluids consisting of polarized particles suspended in a carrier fluid, such as silicone, hydrocarbon based oils, glycol or hydraulic oils. In the case of MR fluids each 3-10 micron diameter ferromagnetic particle has an associated dipole. Upon application of the field, the dipoles align and interparticle forces cause chains to form. Those chains alter the fluid by increasing the viscosity. The shear stress thus becomes field dependent. A key characteristic of the shear stress vs. shear rate or flow curve, is the yield stress. The yield stress is field dependent and defined as the shear stress intercept of the high strain rate asymptote. MR fluids have found applications in dampers, clutches and any other area where a controllable fluid is desired [1–3]. Simulations are a common method of understanding the chain dynamics and how the particle structure changes. Klingenberg has a series of papers on ER fluids [4–6] that covers the basics of many of the interactions in an MR fluid. Ly et. al. present an advanced method for computing the magnetic moment, including multipole effects [7]. Spinks uses simulations to examine chain structure under a variety of fluid flow profiles [8]. Han, Feng and Owen do simulations on the particle dynamics as well as examining the fluid dynamics of the carrier fluid [9]. Mohebi and Jamasbi investigate how particle chains aggregate to form thick columnar structures [10]. Ido examines the effects of non-magnetic particles and spherocylindrical particles [11]. However, these simulations have not been extended to a full MR damper, which has the active fluid in a representative annular channel that is 10mm long, 1mm wide with a radius of 25mm. Representing a full damper would allow for the analysis of the fluid entrance and exit flows, as well as capturing the transient forces and microstructure changes. Considering the rotational symmetry, a representative fluid volume would be 10mm× 1mm× 1mm. With particles of radius a = 3.5μm at a volume fraction of φ = 0.35, a simulation would require approximately 2×107 particles. Existing simulations, using approximately 5000 particles, are too small scale to represent a fluid volume of that size. There are also few good numerical methods for analyzing chain formation at this scale. The standard method for examining the chain formation is looking at particle position time histories are used to show how the chains have formed, as well as the mean length of the chains. While snapshots of particle locations are useful for smaller numbers of particles, they do not scale well to domains having a million particles or more. Similarly, quantities such as mean chain length are dependent on the geometry of the simulation volume, so the values are not scale independent. This paper aims to show the feasibility of simulations at length scales approaching a millimeter and with over one million particles. To achieve this goal, we ran our simulation code on a commercial graphics card, using NVIDIA’s general purpose graphics processing unit programming (GPGPU) environment, CUDA. NVIDIA claims speedups of up to two orders of magnitude over traditional single processor code [12], and we use that ability to run two orders of magnitude more particles than existing simulations [5, 7–11]. To analyze the particle chain formation, we measure the chain length and connectivity, and nondimensionalize the results by comparing it to an ideal case. MAGNETORHEOLOGY MR fluids are magnetic particles suspended in a carrier fluid. Particles are typically made out of iron, cobalt, nickel or an alloy of these. Particle concentration is measured using volume fraction, φ , which is the ratio of particle volume to bulk volume. Under the application of an external magnetic field, H0, particle chains form in the fluid and the fluid thickens. One of the standard descriptions of an MR fluid is as a Bingham plastic, described by τ = τy +ηγ̇ γ̇ ≥ 0, (1) where τ is stress, τy is yield stress, and γ̇ is shear rate. A typical MR device has the fluid flowing through a small millimeter scale gap. There are two typical flow profiles for flow through this gap, a shear mode (Couette flow) and a flow mode, (Pouisieulle flow). In shear mode, which is studied here, the upper surface of the gap is moving, causing a linear flow profile. Assuming a no slip condition at the walls,